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Abstract 

A recently introduced model describing the folding of the triangular lattice is gen¬ 
eralized allowing for defects in the lattice and written as an Ising model with nearest- 
neighbor and plaquette interactions on the honeycomb lattice. Its phase diagram is 
determined in the hexagon approximation of the cluster variation method and the 
crossover from the pure Ising to the pure folding model is investigated, obtaining a 
quite rich structure with several multicritical points. Our results are in very good 
agreement with the available exact ones and extend a previous transfer matrix study. 


PAGS numbers: 05.50.-|-q (Ising problems); 64.60.-i (General studies of phase transitions) 
82.65.Dp (Thermodynamics of surfaces and interfaces). 
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1 Introduction 


Polymerized membranes are two-dimensional generalizations of linear polymers; the large 
variety of possible applications and the extension of one-dimensional statistical properties to 
2D objects justifies recent numerous studies on these systems Differently 

from linear polymers, polymerized membranes are expected to show different long-distance 
behaviours depending on the microscopic characteristics of the system. In particular, rigid 
polymerized membranes are expected to be stable in a flat phase §] which has no analogue 
in polymer systems. 

Models of polymerized membranes consist of networks with fixed connectivity fluctuating 
in an embedding space. The length of the bonds can vary with an energy cost diverging for 
increasing lengths, while the rigidity of the network is described by a bending energy term 
favouring flat configurations. When excluded volume effects are not considered, as it will be 
in the model of this paper, it has been shown 0 B |. B B B that by varying the strength 
K of the bending term, a critical transition arises separating a flat phase at large K from a 
crumpled phase at small K. 

The case of a triangular network embedded in a d-dimensional space where the bonds are 
of hxed length has been hrst considered in |P]. Here the only degrees of freedom correspond 
to the possible ways of folding the network. In P| the problem has been studied in a 2d 
embedding space. In this simplified case the normal vectors to the triangles point either 
up or down in some direction, suggesting a description of the system in terms of Ising spin 
variables. The entropy of this folding problem has been recognized |T^ to be the same 
entropy of the problem of colouring with three colours the bonds of the hexagonal lattice. 
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which has been exactly calculated in . In the description of the model in terms of Ising 


variables a bending rigidity can be easily dehned [H; in [12| a hrst-order 2d folding transition 


between a crumpled and a flat phase has been found at a hnite value of the bending strength 

K. 

In this paper we consider again the 2d folding problem of a triangular network; we apply 
the cluster variation method (CVM) ||I3| , p!5[| to complete the study of [|^ analyzing the 
phase diagram of the model in the whole plane K, h where h is a symmetry breaking held 
distinguishing between up and down triangles. The folding problem can be expressed as a 
vertex model or an Ising problem with a local constraint on the possible spin conhgurations 
T^ . We have also studied the system by progressively relaxing the constraint until the 


standard Ising model is obtained. The relaxing of the constraint corresponds to accept 
folding conhgurations with endpoints in the folding lines, that is cuts between adjacent 
triangles of the network. Therefore intermediate models between the Ising model and the 


folding model of [0, 10, 121 could describe realistic polymerized membranes where defects in 


the connectivity rules are present. 

The cluster variation method is known to be very accurate in describing the phase diagram 


of magnetic systems [1^ . Also in this case of a constrained magnetic system, the CVM gives 
an accurate description of the system. For example, the T = oo entropy calculated with the 
CVM is Ing, g = a/T3/3 ~ 1.2019 which has to be compared with the exact value q = 1.2087 


11|. In the whole plane K, h, in addition to the hrst order transition between the flat and 


the folded disordered phase |^, we have found at negative K a critical transition between 
the disordered phase and a folded phase with staggered antiferromagnetic order. When the 
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constraint is progressively relaxed the topology of the phase diagram is fonnd to evolve in 


qnite a complex way, resembling at negative K the phase diagram of the Ising model for 
metamagnets [jl^, nntil the standard Ising model phase diagram is recovered. 

The paper is organized as follows. In section 2 the model of folding is dehned in terms of 
spin variables and the CVM approximation scheme nsed for stndying the phase diagram is 
briefly described. In section 3 the phase diagrams of the folding system with the constraint 
progressively relaxed are given. In section 4 onr resnlts will be briefly discnssed. 


2 The model and the method 


We briefly describe the model of folding stndied in |p, |^. In a 2d embedding space 

the normal vectors to the triangles of a regnlar network can be represented by Ising spins 
Si = ±1, where i denotes a site of a honeycomb lattice (the dnal of the triangnlar lattice), 
in snch a way that a + (respectively —) spin corresponds to a normal vector pointing np 
(down). Following |^, we consider a bending rigidity K, which measnres the energy cost 
(in nnits of ksT) of a fold between adjacent triangles, and a “magnetic” symmetry-breaking 
held h, conjngate to the normal vectors. Fnrthermore, in order to stndy the crossover from 
the Ising model to the pnre folding model, we shall allow the spins aronnd a hexagon to 
violate the local constraint 


Sj = 0 mod 3 

2 G hexagon 


( 1 ) 


(only 0, 3 or 6 minns spins allowed in a hexagon), which was introdnced in |^. When the 
constraint Eq. |I] is verihed, the model can be expressed as a 11-vertex model |T^. In onr 
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generalized folding model, a violation of the constraint Eq. ^ will cost an energy L > 0, such 
that L = 0 corresponds to the ordinary Ising model, while in the limit L —> cx) one recovers 
the pure folding model. We are thus led to consider the following hamiltonian: 


n 

ksT 


^ ^ ^ SiSj h ^ ^ Si L ^ ^ 6 ('{Sij'ighexagon) i 


( 2 ) 


{ij) i hexagons 

the hrst sum is over nearest neighbors and 5 is a function equal to 1 only when the constraint 
Eq. 2 is satished and zero otherwise; the hamiltonian above will be studied in the hexagon 
approximation of the cluster variation method. 


The CVM, in its modern formulation jl^, [151, is based on the minimization of a free energy 
density functional which is obtained by a truncation of the cluster (cumulant) expansion of 
the corresponding functional appearing in the exact variational formulation of statistical 
mechanics. In the hexagon approximation for the honeycomb lattice the largest clusters 
appearing in the expansion are the hexagons, and the approximate free energy density one 


has to minimize, determined according to [14|, has the form 


/(Pe) = -^it'Tr(siS2P2(si, S 2 )) - hTr(sipi(si)) - pQ{{si}) 

13 1 1 

+ 2Tr(P6lnp6) - ip^^{p2 lnp2) + -Tr(piA Inpi^) + -Tr(pis Inpis) 

+A(Trp6 — 1), (3) 


where Tr stands for trace, Y' denotes summation over the hexagon conhgurations that 
satisfy the constraint Eq. |I|, A is a Lagrange multiplier which ensures the normalization of 
Pe (normalization of the site and pair density matrices follows by dehnition, see below) and 
PiA{B) = Pm(B)(si), P2 = P2(si,S2) and pe = P 6 (si, 52 , S3, S4, S5, se) are the site, pair and 
hexagon density matrices respectively. The spin conhguration of each cluster can be used as 
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a label for the corresponding density matrix, because, as for any classical model, the density 
matrices are diagonal. 

In writing Eq. ^we have introduced two site density matrices, piA and pis, corresponding 
to the two interpenetrating sublattices in which the honeycomb lattice can be divided, which 
will be needed to investigate the case of antiferromagnetic coupling it' < 0. These density 
matrices can be dehned as partial traces of the pair matrix, as follows 

Pia(sa) = y^P2(gA, Sb) 

SB 

Pib{sb) = ^P2(sA,Ss), (4) 

SA 

adopting the convention that the hrst spin in the argument of p 2 always belongs to sublattice 
A. With this assumption, p 2 can be dehned as a (symmetrized, for convenience) partial trace 
of p6 by 

P2{sa,Sb) = ^ [P6isA,SB,s'A,SB,S%s'^) + Pe{s'^,SB,SA,s'j^,S%SB) 

^ c' e" o' c" 

+P6('S(4, s'q, Sa, Sb, s"^, s"q) + P6('S(4, s'^, s"^, Sb, Sa, s"q) 

+P6(sa, s'b, s'^, s"b, Sa, Sb) + Pq{sa, s'b, s'a, s'b, s'a, sb)] , (5) 

where the spins in the argument of pe follow each other counterclockwise in the hexagon, 

and the hrst one is on the A sublattice. In terms of the site and pair density matrices one 

can easily dehne the uniform and staggered order parameters 

^ \ [Tr(spiA(s)) + Tr(spiB(s))] , 

Ms = ^ [Tr(spiA(s)) - Tr(spiB(s))], (6) 

and the nearest neighbor correlation function c = Tr(sASBP 2 (-SA, sb)). 
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With the dehnitions above our free energy can be regarded as a function of pq only and 
taking the derivatives with respect to the generic element of pe we hnd, after some algebraic 
manipulation, the stationarity conditions 


P6(Si,S 2,S3,S4,S5,S6) = exp 


K ^ h ^ 

“ 1 “ r) ^ “ 1 “ TT ^ liGhexagon) 


i=l 


3^^ 


X [P2(S1, S6)P2(S1, S2)P2(S3, S2)p2(s3, S4)p2(s5, S4)p2(s5, 

X [pia{si)pib{s2)pia{s^)pib{sa)pia{s^)pib{s&)]~^^^ , (7) 


where s^ = si and A has to be determined by solving the normalization condition Trpg = 1- 
This set of equations, together with the dehnitions Eqs. | and ^ is known as the natural 


iteration method [13 , because it can be solved by simply iterating the equations above, and 


always converges to a local minimum of the approximate free energy. To hnd the global 
minimum it is therefore enough to start the iteration with diherent sets of initial conditions, 
one for each of the expected phases. 

It can be easily recognized that the elements of pg are not all independent, since when two 
spin conhgurations are related by simmetry (rotation and/or rehection) the corresponding 
elements are degenerate. In Tab. 1 we have listed the 20 independent (up to normalization) 
conhgurations, together with their multiplicities, for the most general case, corresponding 
to the antiferromagnetic phase (see below). In the ferromagnetic phase we have sublattice 
invariance, which implies (using the same symbol for a conhguration and its density) 


^3 — ^45 ^5 — Z &1 ^9 — ^ 10 ) ^11 — ^ 12 , ^15 — ^165 ^17 — ^ 18 ) ^19 — ^ 20 - 

( 8 ) 

In the disordered phase, and for h = 0, we have also invariance under spin inversion, from 
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which one can derive the additional relations 


zi = 2:2, 2:3 = 2^5, = zg, ^9 = ^11, 2:13 = 2:14- ( 9 ) 

Finally, in the limit L —> cx) (pure folding), the densities corresponding to conhgurations 
which violate the constraint Eq. |I| vanish, i.e. Zt = 0, i = 3,4,... 14. 


3 Results 


In this section we present our results for the phase diagram of the model Eq. |^. 

First of all, we consider the case L = 00 , K = h = 0, studied in P, ^ |I^. For 
L = 00 , h = 0 and unbroken spin-flip symmetry the only non-vanishing elements of pe are 
zi = Z 2 , Zi 5 = ziQ, zn = ^18 and 2:19 = Z 2 o- Then the stationarity condition Eq. can be 
solved exactly, yielding 


^19 


+ 3a^ + 6 a + 


a2;i5 = a'^ZiT = a^Zig, 


-1 


a = 


2 — u — \/3 — u — 




u — 1 


u = exp {2K). 


( 10 ) 


For iF = 0 this reduces to 82:1 = 4zi^ = 22:17 = Zig = 4/39, corresponding to an entropy per 

1/3 /2\T2 

site s = Ing with q = 1/13/3 ~ 1.2019, whereas the exact result [pT[| is q = —F ( - ) ~ 

27r V3/ 


1.2087, and to a negative (antiferromagnetic) nearest neighbor correlation c = —1/3, for 
which no exact result is available. The very good agreement (within 0.6 %) of our estimate 
of the entropy with the exact value gives us conhdence in applying the cluster variation 


method to the present model. 









Let us now turn to the analysis of the phase diagrams in the plane K, h at different values 


of L. They are symmetric with respect to the axis h = 0 and it is sufficient to describe them 


at positive values of h. 

In Fig. 1 it is shown the phase diagram of the pure folding problem (L = cx)). At 
sufficiently large values of h the flat phase with M = 1 is always stable. It is remarkable 


that the CVM always yields M exactly equal to 1 in this phase, as conjectured in []^. This 
result is presumably exact, since, as observed in 0, the flat (ferromagnetic) ground state 
has no local excitations. This is no longer true, of course, for hnite values of L. 

The two flat phases with magnetization M = ±1 coexist at h = 0 for K > Kq^{oo) = 
0.1013, to be compared with the value K = 0.11 ± 0.01 of [Q. This hrst order transition 
point can be obtained by requiring that the free energy densities of the ordered and disor¬ 
dered phase take the same value, i.e. by solving the equation /ord(A') = /dis(A')- Here 
/ord(A') = —3iF/2 is the (presumably exact) free energy density of the ordered phase, as in 
, obtained by observing that this phase has vanishing entropy and neglecting the diverg¬ 


ing L contribution, which is common to both phases. On the other hand, by substituting 
the disordered phase density matrix elements Eq. 0 in the expression for the free energy we 
obtain 


/dis(A') 


3 Cr U 1 ,q ^ O \ n 

—K --|- — In Zig 3^19 (o -|- 2cr -|- 2cr) In cr -|- 

2 a+u 2 ^ ’ 


-?,C 


a 

-3/: 

u 

2{a + u) 

2{a -f u) 



In 2, 


( 11 ) 


where C{x) = x\n.x. 

At smaller values of K, the coexistence line at h = 0 separates into two branches where 
the disordered folded phase with M ^ 0 coexists with the flat phases. The intersection of 
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the upper branch with the axis = 0 is at h = 0.18495 ± 0.00005 to be compared with the 


value h = 0.189 found in |]I^. In the disordered phase M vanishes only at h = 0, but this 
might be a consequence of our approximation (see also the discussion in the next section). 

Then, at sufficiently negative value of K, there is a critical transition between the dis¬ 
ordered folded phase and a folded antiferromagnetic phase with staggered order parameter 
Ms 7 ^ 0 and M ^ 0. This transition is represented by the almost vertical broken line of Fig. 
1, which intersects the horizontal axis sX K = —0.284. This completes the phase diagram 


shown in Fig. 6 of [p!2 


In order to understand how this phase diagram evolves while the constraint is relaxed we 
begin by plotting in Fig. 2 (in the L,K plane) the transition line Kq''{L) which separates, 
at h = 0 , the disordered folded phase with M = 0 and the flat phases with M 7 ^ 0 . At 
large values of L the transition is first-order and the curve tends to the asymptotic value 
Kq^{oo) = 0.1013. At small values of L the transition is critical as in the Ising model 
(L=0). At L = 1.359, K = 0.3038 there is a tricritical point where the transition changes its 
behaviour. In our approximation scheme the critical line and the tricritical point are obtained 
as the solution of suitable analytical equations |]^. The derivation of these equations is 


almost straightforward but quite cumbersome, and we have omitted it. 

The topology of the phase diagram in the plane K, h remains the same until the value 
Lb = 1.75 ± 0.05 is reached. For smaller values of L, in the range 1.359 < L < Lb, the 
typical phase diagram is shown in Fig. 3, where the particular value L = 1.6 has been 
chosen. The difference with the case L = 00 is that now the hrst-order transition lines 


between the disordered and the flat phase are interrupted in some range. The two critical 
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points limiting the upper and the lower branch of the interrupted hrst-order line (at positive 
h) have respectively coordinates K = 0.175, h = 0.076 and K = 0.262, h = 0.0009. 

When L still decreases, the branches at smaller \h\ (see the inset of Fig. 3) become 
shorter and shorter until they collapse on the axis h = 0 at the tricritical point of Fig. 2 
(L=1.359). In this situation the phase diagram assumes the topology of Fig. 4 (L = 1.3), 
with the two surviving branches separating phases with different magnetizations. 

By diminishing further the value of L the branches of Fig. 4 become less pronounced 
and at L = 1.1 they are not distinguishable anymore. Therefore, on the transition line 
limiting the antiferromagnetic phase, two tricritical points (in symmetrical positions with 
respect to the K axis) separate the hrst-order behaviour at more negative values of K from 
the critical behaviour close to the K axis. At L = 1.1 the coordinates of the tricritical 
point at h > 0 are K = —0.6198, h = 1.795. These tricritical points move towards larger 
values of |h| when L decreases. It is not clear from our calculations if these points disappear 
at a singular value of L or if they move continuously towards |h.| = cx) when L —> 0. At 
L = 0 the standard phase diagram of the Ising model in a magnetic held is recovered with 
a critical transition always limiting the antiferromagnetic phase |^. This phase diagram is 
shown in Fig. 5. The intersection of the critical line with the axis h = 0 is at iF = —0.6214 
which has to be compared with the exact critical Ising coupling on the hexagonal lattice 


given by \Kc\ = -ln(A/3 -|- 2) ~ 0.6585 |2^. The point K = Kq'^{0) = 0.6214 is the usual 


ferromagnetic Ising critical point. 
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4 Discussion and conclusions 


In this paper we have studied the phase diagram of a constrained Ising spin model describing 
a network of equilateral triangles embedded in a 2d space. The constraint is due to the fact 
that during the folding the local connectivity properties of the network have to be preserved. 
This study has been done by applying the cluster variation method to the hamiltonian defined 
in Eq. In Eq. § the term proportional to L favours the spin conhgurations verifying the 
constraint which are the only surviving in the limit L ^ oo. In order to study networks with 
defects consisting of endpoints in the folding lines, we have also studied the model Eq. || at 
finite values of L by progressively relaxing the constraint until to consider the case L = 0 
corresponding to the usual Ising model. 

Our results can be summarized in this way. In the case of the pure folding problem 
we first observe the good agreement between our evaluation of the entropy s = Ing with 


q ~ 1.2019 at = 0, h = 0 and the exact result q = 1.2087 of |^. We also predict for 
the nearest neighbor correlation at K = 0, h = 0 the value c = —1/3. Results concerning 
other regions of the phase diagram in the plane K, h are also in good agreement with results 
obtained in |]^. We have completed the results of [|T^ by studying the phase diagram of the 
model Eq. also for negative values of the bending rigidity K. A critical transition between 
the disordered folded phase and a folded phase with antiferromagnetic order has been found. 
The critical value of K for this transition at h = 0 is A' = —0.284. We conclude the discussion 


of the pure folding case making some comments about the hypothesis advanced in [12] that 


M always equals 1 in the ordered phase and 0 in the disordered phase. Our results conhrm 


that M = 1 in the ordered phase, a result which is almost certainly exact, but yield M = 0 
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in the disordered phase only for h = 0, while close to the transition with the flat phase 
we have M ~ 10“^. An order parameter of the same order of magnitude can be obtained 
by extrapolating with standard methods (Shanks transform, alternating e-algorithm and 
Fade approximants [^, the results from the transfer matrix method proposed in [VI 


obtained with strip widths in the range 2 to 8. It would be very interesting to see how the 
order parameter varies when one considers larger clusters for the CVM or larger strips in 
the transfer matrix method, but this is beyond the purpose of the present paper. In our 
opinion, however, this issue cannot be settled numerically in a dehnitive way, and a rigorous 
argument would be welcome. 

At hnite L, the evolution of the phase diagram from L —oo to L = 0 is shown in Figs. 
1,3-5. It is curious to observe the analogies between the antiferromagnetic region of these 
phase diagrams and the phase diagrams appearing in the study of the Ising spin model for 
metamagnets [^. The hamiltonian of the Ising spin model for metamagnets includes an 
antiferromagnetic coupling J between nearest neighboring sites, a ferromagnetic coupling 
J' between next-to-the-nearest neighboring sites and an external magnetic held. When the 
ratio e = J j J' is in the range 0 < e < 3/5, the phase diagram of the metamagnet is similar 
for example to that of Fig. 4, but with the hrst-order branches inside the antiferromagnetic 
phase. The branches terminate at critical points where two antiferromagnetic phases with 


different net magnetizations cease to coexist |17]. These analogies can be understood by 


observing that the constraint energy proportional to L includes an effective next-to-the- 
nearest neighboring site interaction in the hexagonal lattice. 


In conclusions, we observe that in this paper the CVM has been applied to a vertex 


13 

















model giving results in excellent agreement with numerical and exact previously known 


results. This is encouraging for applying the CVM to the study of other vertex models or 


colouring problems. In particular, in |23] an extension of the model of P, 10, 121 which 


describes a possible embedding of a triangular network in a 3d space has been formulated as 


a 98-vertex model. We think that the phase diagram of the model of could be efficiently 
studied by applying the CVM, and work is in progress along these lines. 
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Figure Captions 


Fig. 1 : Phase diagram of the pure folding problem (L = oo). The solid and the broken 
lines are respectively coexistence and critical lines. 

Fig. 2 : The transition line separating at h = 0 the disordered folded phase (M = 0) from 
the flat phases (M = ±1). The solid and the broken lines are respectively first-order and 
critical lines. The tricritical point is marked by a full circle. 

Fig. 3 : Phase diagram of the folding model at L = 1.6. The solid and the broken lines 

are respectively coexistence and critical lines. The full circles represent critical points (their 

coordinates are given in Section 3). In the inset it is reported the magnihcation of the region 
of the phase diagram in the empty square. 

Fig. 4 : Phase diagram of the folding model at L = 1.3. The solid and the broken lines 

are respectively coexistence and critical lines. The full circles represent critical points; their 

coordinates are: K = —0.1382, h = 0.66; K = —0.1382, h = —0.66; K = 0.31377, h = 0. 

Fig. 5 : Phase diagram of the pure Ising model on a honeycomb lattice (L = 0). The solid 
and the broken lines are respectively coexistence and critical lines. The full circle represents 
the usual ferromagnetic Ising critical point. 



Table Caption 


Table 1: Independent hexagon configurations. 



Table 1 


Configuration 

Si 

-S2 

S 3 

S 4 

S 5 

S6 

Multiplicity 

2:1 

+ 

+ 

+ 

+ 

+ 

+ 

1 

2:2 

— 

— 

— 

— 

— 

— 

1 

^3 

— 

+ 

+ 

+ 

+ 

+ 

3 

^4 

+ 

— 

+ 

+ 

+ 

+ 

3 

^5 

+ 

— 

— 

— 

— 

— 

3 

2^6 

— 

+ 

— 

— 

— 

— 

3 

2:7 

— 

— 

+ 

+ 

+ 

+ 

6 

2^8 

+ 

+ 

— 

— 

— 

— 

6 

2:9 

— 

+ 

— 

+ 

+ 

+ 

3 

ZlO 

+ 

— 

+ 

— 

+ 

+ 

3 

Zll 

+ 

— 

+ 

— 

— 

— 

3 

Zl2 

— 

+ 

— 

+ 

— 

— 

3 

Zl-i 

— 

+ 

+ 

— 

+ 

+ 

3 

Zli 

+ 

— 

— 

+ 

— 

— 

3 

Zlh 

— 

— 

— 

+ 

+ 

+ 

3 

ZlQ 

+ 

+ 

+ 

— 

— 

— 

3 

zn 

— 

— 

+ 

— 

+ 

+ 

6 

Zl% 

+ 

+ 

— 

+ 

— 

— 

6 

Zl 9 

— 

+ 

— 

+ 

— 

+ 

1 

Z20 

+ 

— 

+ 

— 

+ 

— 

1 
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